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ABSTRACT 


Propagation  of'  electromagnetic  energy  through  the  atmosphere  is  a  diHicult  task 
because  of  temperature  fluctuations  and  index-of-rcfraction  inhomogeneities  which  de¬ 
grade  the  beam's  coherence.  Understanding  this  phenomena  is  of  practical  importance 
for  optical  systems. 

This  thesis  presents  an  analytical  numerical  technique  which  simulates  the  effects 
of  atmospheric  turbulence.  The  extended  Huygens- Fresnel  principle  was  used  to  simu¬ 
late  wave  propagation  in  a  two-dimensional  randomly  vary  ing  medium,  which  is  repres¬ 
ented  by  thin,  filtered,  Gaussian  phase  screens.  The  wave  optics  code  implements  both 
Fresnel  and  Fraunhofer  propagation,  by  employing  the  fast  Fourier  transform  (FFT) 
algorithm.  The  analytical  spatial  coherence  length,  ,  and  noimalizcd  intensity  vari¬ 
ance,  (7-//-,  of  the  perturbed  electric  field,  were  examined.  Results  support  the  concept 
of  intensity  saturation  for  weak  scattering  cases,  however.  dilTerences  in  the  values  of  the 
theoretical  and  analytical  spatial  coherence  lengths,  occurred. 
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I.  INTRODUCTION 


Atmospheric  turbulence  degrades  the  coherence  of  electromagnetic  energy  propa¬ 
gating  through  the  atmosphere.  The  refractive-index  inhomogencitics  associated  with 
the  turbulent  atmosphere  induce  phase  and  intensity  perturbations  across  the  wavefront. 
Understanding  the  propagation  and  scattering  of  optical  waves  in  random  media,  is  es¬ 
sential  for  atmospheric  laser  beam  propagation  and  imaging  systems. 

This  thesis  models  the  propagation  through  a  random  media  by  means  of  the  ex¬ 
tended  Huygens-Fresnel  principle.  A  two-dimensional,  thin,  Gaussian  phase  screen  re¬ 
presents  the  randomly  varying  medium.  This  wave  optics  propagation  code  employs  the 
fast  Fourier  transform  (FFT)  algorithm  in  both  Fresnel  and  Fraunhofer  forms  of  prop¬ 
agation.  The  Fresnel  region  incorporates  two  forms  of  the  diffraction  integral,  the 
transfer  function  form  for  "near  field"  distances,  and  the  convolution  form  for  "far  field" 
regions. 

Previously,  numerically  intensive  simulations  of  this  type  required  large,  super¬ 
computer  systems.  However,  these  computations  were  performed  on  a  compact, 
desktop  computer  system. 

The  model's  accuracy  was  assessed  by  comparing  computer  values  of  the  spatial 
coherence  length,  .  with  values  obtained  from  the  analytical  mutual  coherence  func¬ 
tion  ( MCF).  .Additionally,  the  concept  of  intensity  variance  saturation  was  examined  for 
a  single  phase  screen  in  the  Fraunhofer  approximation.  1  he  normalized  intensity  vari¬ 
ance.  appr: aches  a  saturation  value  asymptotically  for  increasing  values  of  the 

index-of-refraction  structure  parameter,  Q  .  Theoretical  calculations  suggest  saturation 
for  multiple  scattering,  however,  they  say  nothing  about  a  single  phase  screen  realiza¬ 
tion.  This  thesis  provides  results  which  support  saturation  clfects  for  single  scattering 
cases. 
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II.  BACKGROUND 


A.  STATISTICAL  DESCRIPTION  OF  TURBULENCE 

1.  Random  Variables 

Maintaining  the  coherence  of  an  electromagnetic  wave  propagating  through  the 
atmosphere,  requires  an  understanding  of  the  effects  imposed  on  the  wave  by  the  tur¬ 
bulent  medium.  The  fundamental  characteristic  of  atmospheric  turbulence  is  its  ran¬ 
domness.  which  must  be  described  statistically.  In  addition  to  statistical  quantities, 
further  assumptions  must  be  made  of  atmospheric  turbulence.  These  include  the  con¬ 
cepts  of  stationarity.  homogeneity  and  isotropy.  Stationaritv  implies  tliat  random 
processes  are  time  independent.  Homogeneity  assumes  invariance  under  a  Galilean 
transformation  of  coordinates,  while  isotropic  variables  are  invariant  with  respect  to 
coordinate  rotations.  (Ref.  1]  Mathematically,  these  two  assumptions  imply  that  the 
statistics  at  two  points  and  depend  only  on  the  dilTcrence.  r,,  =  if,  —  f,!  . 

2.  Local  Homogeneity  and  Isotropy 

In  general,  atmospheric  random  variables  do  not  obey  the  assumptions  of 
stationarity.  homogeneity,  and  isotropy.  However.  Tatarski  (Ref.  2)  introduces  the 
concept  of  "local"  iiomogeneous  and  isotropic  random  variables.  Tins  concept  requires 
homogeneity  and  isotropy  within  a  localized  region  of  size  Z-,,.  the  outer  scale  length. 
Furthermore,  th.e  diiTiculties  associated  with  nonstatiotiars  random  variables,  are  re¬ 
moved  by  considering  random  lields  with  stationary  first  increments  [Ref.  2).  Tatarski  s 
meiliod  applies  to  a  nonstationary  random  function  whose  mean  varies  slow  ly  with  time, 
by  considering  the  dilfcrencc  of  the  {'unction  at  two  difierent  locations.  1  he  slow  func¬ 
tional  changes  do  not  alTect  the  \  aluc  of  this  dilfcrencc. 

3.  Structure  Functions 

Tatarski  introduces  the  structure  function 

D/T2  -  fi)  =  <1/1'':)  (1) 

a  tensor  that  is  the  difference  between  two  quantities.  Some  important  aspects  of  the 
structure  function  are;  that  its  general  form  is  valid  for  any  variable,  and  <  >  denotes 
an  ensemble  average  taken  over  all  possible  point  pairs  f,  .  T,  Assuming  homogeneity 
and  isotropy,  the  vector  dependence  reduces  to  the  magnitude  ri2=|r2-f|i  .  and  the 
structure  function  becomes 


/)/(/■)- <i/(''2) -yi'i)r>. 


Kolmog  .ov  showed  through  dimensional  analysis,  a  simple  power  law  dependence  of 
equation  (2).  over  a  limited  interval  called  the  inertial  subrange,  as 

Dj(r)  =  cjy^'\  0) 

Tatarski  introduces  the  concept  of  "passive  additives",  which  are  quantities  independent 
of  position  in  the  vector  field,  and  do  not  directly  influence  the  dynamics  of  the  turbulent 
medium.  Temperature  is  a  passive  additive  and  has  a  structure  function  of  the  form 

D^{r)  =  (4) 

As  long  as  r  remains  within  the  inertial  subrange,  temperature  is  approximated  as  a 
passive  additive,  and  equation  (4)  is  valid.  Likewise  the  index-of-refraction  is  a  passive 
additive  with  a  structure  function 

/),,(;•)  =  (5) 

Since  tlie  index  of  refraction  depends  on  the  density  of  the  atmosphere.  D,  and  D,  are 
related  by 


4.  Coiarianee.  Pmser  Spectral  Densities 

In  addition  to  the  structure  function,  other  characteristics  of  random  processes 
include  the  covariance  lor  correlation),  and  Power  Spectral  Densities.  It  is  the  interre- 
lation'>hip  of  theve  three  quantities  wiiich  provide  a  useful  method  for  analysing  random 
processes.  The  covariance  between  two  random  variables  S  and  T  can  be  expressed  as 

/isr=  -  <7'(r.)>]l.'>(r,)  -  <S(r,)>]>.  (7) 

However,  more  frequently  it  is  the  autocovariance  function 

^/T=  <7(/|)>linr:)  -  <T(r.}>]>.  (S) 

which  is  needed.  Luitlicrmcire.  if  T  is  homogeneous  and  <T>  —  0  is  assumed,  then 
equation  iS)  simplifies  lo 


[^)] 

Combining  this  relation  and  equation  (2).  gives  an  expression  for  the  relationship  be¬ 
tween  the  structure  function  and  covariance  function  as 

Djir)  =  2\B-j-j{^)  —  Bjjir)].  (10) 

In  one  dimension  the  covariance  function  and  the  power  spectral  densitx'  are  transform 
pairs  given  by 


and 


n'(Kj  = 


(in 


/>(  n 


f'+y. 

e‘’^''ir{KUfK. 


(12) 


I'sing  the  fact  that  B(r!  is  an  even  function  and  substituting  equation  (12)  into  equation 
( I'O.  D,n  becontes. 


/)(<•)==  2j  (1  -  cos(Kr)l(r(KWK-.  (13) 

Tatarski  devei,  ps  an  expression  I'or  the  one-dimensional  Kolmogorov  spectral  density, 
giver,  by 

/r(K)  =  0.1224f;>"^'-\  (14) 

Discussion  to  this  point  has  been  of  one-dimensional  random  processes,  how¬ 
ever  these  concepts  arc  applicable  to  three-dimensional  cases.  Analogous  to  equation 
(11).  Tatarski  defines  the  three-dimensional  power  spectral  density  as 


4 


•+OC 


/'+OC  /*+oo 


(l?iK)  = 


— oo  — oc  — oo 


and  similarly,  tlie  correlation  function  is 


B{r) 


•+00 


(16) 


Using  the  relation 


(D(k') 


_1  dJJ'U) 
2~k  c{k 


the  three-dimensional  Kolmogorov  power  spectral  density  becomes 

(D(k)  =  0.033Qk“'''1 


(H) 


(IS) 


B.  EM  PROPAGATION  THROUGH  TURBULENCE 
1,  Wave  Equation 

Based  upon  Tatarski.  ClilTord  [Ref.  3)  develops  theoretical  results  oflinc-of-siglit 
propagation  through  the  atmosphere,  directly  from  Maxwell's  equations.  Assuming  zero 
conductivity,  and  unit  magnetic  permeability  in  the  atmosphere,  as  well  as,  a  sinusoidal 
time  dependent  electromagnetic  field.  Maxwell's  four  equations,  in  Gaussian  units,  take 
the  form 


V  .  //  =  0, 

(19) 

<1 

X 

II 

(20) 

Vx//  = 

(21) 

V  .  (,E£)  =  0. 

(22) 

Taking  the  curl  of  equation  (20),  and  substituting  it  into  equation  (21).  gives 

-  T-  V(T .  T\  =  /<  'n ( 23 ) 
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Rewriting  equation  (22)  in  the  form 

£.V«^  +  //^V.£  =  0,  (2-lj 

and  substituting  it  into  equation  (23).  yields  the  vector  form  of  the  wave  equation 

V‘£  + + 2V(£.  Vlogn)  =  0.  (25) 

The  third  term  of  equation  (25)  describes  the  change  in  polarization  of  a  propagating 
electromagnetic  wave.  This  term  is  negligible  as  long  as  the  wavelength  is  small  com¬ 
pared  to  the  refractive  inhomogeneities.  Thus  equation  (25)  reduces  to 

v"E  -f  k^n^E  =  0.  (26) 

Equation  (26)  is  the  vector  form  of  the  wave  equation  describing  propagation  through 
the  turbulent  atmosphere.  The  difficulty  in  solving  this  equation  lies  in  the  second  term 
containing  the  random  variable  n.  Various  methods  are  available  for  obtaining  solutions 
to  equation  (26),  each  of  which  relics  on  several  critical  approximations.  Strohbelm 
[Ref  4]  lists  these  approximations  as: 

1.  Negligible  depolarization  elTects. 

2.  Negligible  back-scattering. 

3.  Use  of  the  parabolic  approximation  to  the  wave  equation. 

4.  Turbulence  is  uncorrelated  in  the  diiection  of  propagation. 

2.  The  Method  of  Small  Perturhations-Boni  Approximation 

Both  Tatarski  [Ref  2]  and  ClilTord  [Ref  ?].  solve  tlic  wave  equation  in  a  turbu¬ 
lent  atmosphere  using  the  method  of  small  perturbations,  which  is  equivalent  to  the 
Born  approximation.  This  method  expands  the  electric  field  into  a  series  of  decreasing 
amplitudes,  and  the  refractive  index  into  a  power  scries  in  the  form 


£  —  £c  +  £]  "1"  — 

(27 

«  =  1  4-  -f- . 

(28 

Substituting  these  into  equation  (26)  and  equating  same  order  terms,  results  in  two 
equations 

V‘£,j -f  A*£,  =  U.  (26) 
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+  A-^£i  +  2A-«,Eu  =  0. 


(50) 


where  terms  of  order  /jJ  and  higher  are  ignored.  Assuming,  as  Tatarski  does,  that  tlie 
unperturbed  field  is  a  plane  wave  propagating  in  the  z-direction  represented  as 
fo  =  exp[ikz]  ,  equation  (30)  becomes 

+  =  (31) 

an  inhomogeneous  partial  differential  equation  with  constant  coefficients.  Its  solution 
is  the  convolution  of  a  plane  wave  Green's  function  with  the  source  term,  or  inhomo¬ 
geneous  term,  given  by 


e\p(/A;ir  —  I'i) 

d-i  ' - - — — - [2/v-/i|()-')  e.vpf/A-r')]. 


|r  —  / 


3.  The  Metliod  of  Smooth  Perturbations-Rylov  Appro.ximation 

In  addition  to  the  Born  appro.ximation,  Tatarski  develops  the  Rytov  approxi¬ 
mation  ivhich  assumes  a  solution  to  equation  (3(>)  of  the  form 


E  =  expl'l')  =  exp(.\  +  iS) 


or  simply 


r.  =  ,1  exp(/5) 


34) 


where  .A  is  the  amplitude  given  by  A  -  expX  .  Applying  equation  (33)  to  equation  (30) 
and  dividing  bv  vicld<;  the  Rviov  solution  eiven  bv. 


V‘E  4-  k^n^(r)  =  V  log  E -f  (V  log  EV  4-  k'n"{r). 


(35) 


Tatarski  farther  shows  that  both  methods  of  approximation  arc  equivalent. 

C.  HUVGENS-FRESNEL  THEORY 

As  we  have  seen,  the  theories  of  Tatarski  and  Clilford  use  the  dilTerential  equation 
approach  to  solve  the  problem  of  propagation  through  a  turbulent  medium.  However, 
Lutomirski  and  Yura  [Ref  5).  approach  this  problem  in  terms  of  integral  equations 
which  use  an  extended  lluygens-Frcsnel  theory.  This  technique  is  equivalent  to  a  dif¬ 
ferential  equation  approach,  but  it  is  easier  to  integrate  and  simulate  using  m  tech- 


niques.  Lutomirski  and  Yura,  develop  an  extended  Huygens-Fresnel  theory  by 
introducing  a  random  phase  term  for  turbulence  in  the  Huygens-Fresnel  integral  which 
is  developed  in  standard  optics  texts  like  [Ref.  6].  This  additional  phase  perturbation 
takes  the  form  of  the  Rytov  approximation,  .  The  extended  Huygens-Fresnel  integral 
is, 


-ik 

2n 


e(/A:|F-  r'l) 


(36) 


In  the  geometrical  optics  limit,  Fermat's  Principle  is,  '¥-^k f  /2,(z)c/z.  With  a  power  series 
expansion  of  c'*',  equation  (36)  reduces  to  the  Green's  function  solution  of  Tatarski  and 
Clifford  given  in  equation  (32). 

Recently,  Martin  and  Flatte  [Ref.  7]  presented  an  atmospheric  turbulence  algorithm 
which  uses  the  dilTcrential  equation  approach  that  has  as  a  solution,  the  extended 
Huygens-Fresnel  integral.  A  filtered  random  Gaussian  phase  screen  introduces  the 
phase  perturbations  while  the  algorithm's  path  integral  method,  incorporates  a  multi¬ 
screen  transfer  function  form  of  Fresnel  propagation. 

D.  MUTUAL  COHERENCE  FUNCTION 

The  effects  of  atmospheric  turbulence  can  be  expressed  in  terms  of  two  functions, 
the  Mutual  Coherence  Function  (MCF)  and  the  Modulation  Transfer  I'unction  (MTF). 
Lutomirski  and  ^'ura  [Ref.  5]  derive  the  first  concept  by  considering  the  average  intensity 
</(/■»  =  ,  of  equation  (36).  What  results  is  an  at  crage  intcn.sity  which  is 

a  product  of  the  autocovariance  of  the  aperture  and  the  atmospheric  MCF.  The  atmo¬ 
spheric  .MCF  term  has  the  Rytov  form  <  expiT" -h 'l-"')>  where  the  T"  refers  to  the 
complex  phase  factor  at  the  i-j  coo.’-dinatc  and  'F  '  '  corresponds  to  the  r-  coordinate.  This 
term  is  log-normally  distributed,  as  Jong  as,  M'  is  composed  of  Gaussian  variables.  Using 
this  fact  and  results  in  Fried's  work  [Ref  S  ],  the  atmospheric  MCF  was  written  in  terms 
of  the  wave  structure  function  D(/7  ).  given  by 


-exp[.^],  ,37) 

where  Dip)  =  Dyip)  +  D^ip)-  Lutomirski  and  "^'ura  apply  the  structure  function  for  a 
plane  wave. 
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Cl{z)dz,  (3S) 

to  equation  (37)  which  can  be  written  as, 

-  MCFU,) .  exp[-(  ^  )’■"].  (39) 

where  Po  ,  the  lateral  cohernece  length,  given  by 


Z)(p)  =  2.91  A 


Po 


AGk 


Q(z)./r] 


-3/5 


where  k  =  and  C-  is  the  index  of  refraction  structure  parameter  along  the  optica! 
path  length  L.  represents  the  distance  where  the  spatial  coherence  of  the  wave  drops  to 
c"'  point  of  the  MCF.  [Ref  8) 


E.  MODULATION  TRANSFER  FUNCTION 

Lutomirski  and  Yura's  concept  of  MCF  is  closely  related  to  the  MTF  of  Fried's 
[Ref.  S].  The  MCF  and  MTF  arc  in  fact  the  same  function  but  expressed  in  terms  of 
dilfercnt  variables.  The  MCF  is  measured  in  the  coordinates  of  the  propagation  field 
and  has  dimensions  of  distance,  while  the  MTF,  is  measured  in  the  image  plane  and  has 
the  dimensions  of  spatial  frequency.  Both  are  equivalent  under  a  transformation  be¬ 
tween  the  two  planes  by  letting  p-*/./?/ where  R  is  the  focal  length  of  the  optics  and  f  is 
the  spatial  frequency.  This  transformation  is  valid  under  the  Wiencr-Khintchine  theo¬ 
rem  since  the  lens  Fourier  transforms  the  incident  electric  field  at  the  aperture  to  the 
image  plane.  Fried's  expression  for  the  atmospheric  long  term  .MTF  is 


.\!TFif)  =  exp! 


MV/3 

''0  j 


(41) 


where  =  2. 1  pj. 

In  calculating  the  .MTF,  two  distinct  cases  exist.  These  are  a  short  term  and  long 
term  MTF.  d  he  short  term  exposure  describes  the  evaluation  of  the  wavefront  in  sufFi- 
ciently  short  time  intervals  such  that  the  turbulence  appears  frozen.  The  long  term 
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MTF,  is  a  single  long  time  integrated  exposure,  taking  into  account  every  turbulence 
configuration.  This  thesis  focuses  on  the  method  prescribed  by  the  short  term  .M'ff. 
[Ref  S] 

The  analytic  distinction  between  the  two  cases  lies  in  the  manner  in  which  the 
wavefront  distortions  are  handled.  Specifically,  the  distortion  attributed  to  a  random  tilt 
of  the  wavefront.  Tilt  results  from  the  varying  phase  fluctuations  across  the  aperture 
which  accumulate  along  the  optical  propagation  path.  These  fluctuations  produce  im¬ 
age  motion  in  the  focal  plane  of  the  receiver.  For  a  very  short  exposure,  the  tilt  is  ex¬ 
tracted,  by  fitting  a  mean  square  two-dimensional  flat  plane  to  the  electric  field  and 
rotating  it  through  an  angle  so  that  the  mean  wavefront  is  normal  to  the  direction  of 
propagation.  This  introduces  a  phase  shift,  resulting  in  the  displacement  of  each  curve 
about  the  optical  axis.  Fried's  [Ref  S]  development  of  this  theory  suggests  a  higher 
MTF  at  all  spatial  frequencies  for  the  short  term  MTF.  The  short-exposure  M  IT  for 
near-field  and  far-field  cases  is  respectively  given  by 


<'oi{)>n/=  exp 
-roOO  exp 


r  ar/i 

5  3 

- 

r 

;  '2 

-3.44 

''o 

1  - 

1 

L  . 

-3.44 


^■Rf 


5/3 


/.Rf 

D 


1/3 


(-12) 


where  Tq  is  the  MTF  of  a  diffraction-limited  lens,  and  the  exponential  term  corresponds 
to  the  atmospheric  MTF.  These  equations  predict  a  near  field  short  term  M  f F  that 
start  at  one.  declines  to  a  minimum,  then  increases  to  unity,  at  the  optical  cut  olf  fre¬ 
quency.  Figure  I  illustrates  this  phenomena. 

F.  DIFFR.ACTION  INTEGRAL 

This  section  presents  the  analytical  work  concerned  with  the  development  of  my 
propagation  algorithms.  It  begins  by  illustrating  the  approximations  that  were  used  to 
manipulate  the  diflraction  integral  and  then  proceeds  with  the  analysis  required  to 
transform  the  solution  of  the  I  lelmbioltz  equation  into  two  forms  of  the  Huygcns-Fresnel 
principle.  The  two  forms  are,  first,  a  convolution  form  suited  for  long  distance  propa¬ 
gation  and  second,  as  a  transfer  function  form  for  short  distance  propagation.  Some 
initial  assumptions  made  by  Roberts  (Ref  9J  include  the  following: 

1.  Light  propagates  in  the  k  direction. 

2.  The  wave  amplitude  is  known  in  the  xy-piane  at  z  =  0  . 
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Mutual  Coherence  Function 


Figure  1.  Short  Term  Mutual  Coherence  Function  for  an  8  x  8  Subaperture. 
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3.  Polarization  efTects  are  negligible. 

4.  The  electric  field  amplitude  is  a  scalar  function  V(r,z). 

Following  Roberts,  the  analysis  begins  with  the  Huygens-Fresnel  integral  for  the 
propagation  of  light  waves  given  as 

y(F,2)  =  ~f  f^V(p,  0)  exp[^  v'z'  +  1^  -  Pi'  ]• 

where  /  is  the  wavelength  of  light  and  p  is  a  vector  in  the  aperture  plane,  Fisa  vector 
in  the  image  plane,  and  z  is  the  propagation  direction.  From  the  Fresnel  approximation 
where,  r  <  <  z  and  p  <  <  z,  factor  from  the  square  root  and  expand  by  a  binomial 
expansion.  Equation  (43)  becomes 


Since  the  exponential  phase  factor  outside  of  the  integral  does  not  affect  intensity 
measurements,  equation  (44)  is 

F(r.  r):s;  — ^  JoT'fXp.  0)  exp|^  |F  -  pi' j-  (-1-) 

Expanding  the  quadratic  term  gives, 

1  (F,  z)  =  ^jdpf  ’(p,  0)  expl”  [r^  -  2F  •  p  /  p^]l, 

^  ^  r  -  -  1  (■^^’) 

=  -^cxp|^-^r^JJt/pr(p,0)expj^-g-p^J  expj^  -27r/j^-^^ 

This  is  the  convolution  form  of  the  Fresnel  integral,  which  is  equivalent  to  the 
Fraunhofer  integral  except  for  the  quadratic  phase  factor  in  the  integral. 

To  obtain  the  transfer  function  form,  the  analysis  begins  with  denoting  as  the 
Fourier  transform  of  \’(r,z)  given  by 
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*’(/'.-)  =  I  J^exp(-2n-//.r)r(r,  r),  (47) 

where  l  '{r,  z)  is  given  in  equation  (45).  Interchanging  the  order  of  integration  yields 

I  '{f,  z)  =  -rf  JdT)  r(p,  0)1^/;"  exp(  -2nif.r)  expj^  “fr 

Next,  a  change  of  variables  is  made  where 

r'  =  r-  p.  (49) 

Substituting  this  into  equation  (4S)  gives  the  following  equation 

I  'if,  z)  =  -jy  Jt/p  f  (p.  0)Jt/r'  exp[  —2nif  •  [r  +  p)]  expF  r'^1. 

-■  c  -  r  .  r  •  -1 

=  — f  J  >  '(p,  0)  exp(  -2?://' .  p)J  dr'  exp(  -2-if .  P)  expj^  “fp  J' 

The  r '  integral  is  replaced  with  its  Gaussian  transform  pair. 

u:c\p[~i-).z/).  (51) 

giving 

I  ■(/’.  z)  =  exp(  -in/.zf)^ (f,  J  '{p.  0)  cxpl  -2r/f .  p ).  (52 ) 

1  he  inverse  Fourier  transform  of  equation  (52)  yields 

I  V".  t:)  =  |t/7exp(27r// .  r)  exp(  -in/.zf)jifp\'{p,  0)  exp(  -2nif .  p).  (53) 

This  expression  is  the  transfer  function  form  of  the  dilfraction  integral,  which  is  equiv¬ 
alent  to  the  solution  of  the  differential  equation  approach  used  by  Martin  and  ITatte 
[Ref  7J.  Reviewing  the  form  of  equations  (46)  and  (53).  the  need  for  two  equations  de¬ 
scribing  different  propagations,  is  obvious.  In  one  instance,  the  propagation  distance  z 
enters  in  the  denominator  of  the  exponential  term.  It  is  at  long  distances  that  the  ex¬ 
ponential  term  varies  slowly.  On  the  other  hand,  equation  (53)  is  suited  for  short 
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propagation  distances,  as  2  enters  into  the  numerator  of  the  exponential  for  slow  vari¬ 
ations  in  this  case. 
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III.  NUMERICAL  SIMULATION  .MODEL 

This  numerical  simulation,  modeled  wave  optics  propagation  of  an  electromagnetic 
wave  through  a  random  medium,  represented  by  a  two-  dimensional  Gaussian  phase 
screens.  Techniques  in  this  model  required  certain  "tools"  and  their  testing.  Ihese 
"tools"  included  a  need  for  a  reliable  random  number  generator,  used  in  generating  the 
random  phase  screens,  and  an  efficient  two-dimensional  FFT,  used  in  approximating  the 
diffraction  integrals.  But  first,  a  discussion  of  the  experimental  arrangement  is  needed. 

A.  EXPERIMENTAL  ARRANGEMENT 

Due  to  the  extensive  numerical  calculations  in  this  simulation,  a  Compaq  deskpro 
S03S6-20  computer  was  used.  It  features  a  64  megabyte  hard  drive  and  16  megabytes 
of  memorx'.  In  addition  to  the  20-MMz  803S7  coprocessor,  a  Weitek  1167  math 
coprocessor  was  added  to  enhance  execution  speed.  The  32  bit  Fortran-3S6  compiler 
was  from  Silicon  Valley  Software  (SVS)  and  uses  Phar  Lap  Software  to  extend  the  op¬ 
erating  .system  beyond  one  megabyte.  The  math  and  graphics  packages  were  produced 
by  Scitcch  Scientific.  The  Compaq  has  a  640  x  4S0  pixel  VG.A  grapliics  monitor.  Both 
the  HP  LaserJet  Scries  II  and  Panasonic  K.\-P1092i  multi-mode  printers  were  used  in 
this  arrangement. 

B.  CO.MPUTER  PRELIMINARIES 
1.  Random  number  generator 

The  main  purpose  of  a  computer  simulation  is  to  approximate  natural  phe¬ 
nomena.  To  make  things  realistic,  random  number  sequences  were  used  to  iinroducc 
stochastic  variations.  One  might  ask,  wlutt  minimal  criterka  should  a  particular  random 
number  generator  satisfy?  Certainly  the  most  important  criterion  is  that  the  sequence 
of  numbers  is  sufficiently  random.  Other  criteria  are  uniformity,  reproducibility,  mini¬ 
mum  memorx',  fast,  non-repeating,  and  statistically  independent. 

The  more  difficult  characteristic  to  satisfy  is  statistical  independence.  Thus  a 
series  of  tests  were  needed  to  provide  a  quantitative  measure  of  the  generator's  per¬ 
formance.  There  are  two  kinds  of  statistical  tests;  empirical  tests  and  theoretical  tests. 
Empirical  tests  focus  on  how  the  computer  manipulates  groups  of  numbers  from  the 
sequence  and  evaluates  certain  statistical  quantities.  Perhaps  the  best  known  of  all  sta- 
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tistical  tests  are  the  "Chi-Square"  tests.  Theoretical  tests,  on  the  other  luind.  establish 
characteristics  of  a  sequence  using  methods  based  on  recurrence  rules.  [Ref.  10] 

Tor  the  purpose  of  this  simulation,  only  empirical  tests  were  applied  to  the 
SVS-Fortran  random  number  generator  called  RanH j  .  The  following  five  tests  were 
considered; 

1.  Frequency  Test.  This  test  determines  whether  or  not  the  sequence  of  numbers  are 
uniformly  distributed  as  L’(O.l). 

2.  Serial  Test.  This  test  is  an  extension  of  the  frequency  test  to  two  dimensions  or 
matrix  form. 

3.  Fagged-Product  Test.  This  test  checks  for  correlations  between  successive  numbers 
over  a  given  lag  period. 

4.  Run  Tests: 

a.  Runs  up  and  down.  This  tests  for  long  increasing  and  decreasing  sequences  of 
numbers. 

b.  Runs  above  and  below  the  mean.  This  tests  for  long  sequences  with  values 
consistently  abuse  or  below  the  mean. 

The  results  of  the  five  tests  arc  provided  in  Appendix  A.  Of  these  five  tests,  the 
Fagged-Product  and  Runs  Tests  arc  the  most  critical  when  simulating  atmosplieric  tur¬ 
bulence.  Tiiesc  three  tests  determine  wlicther  or  not  correlation  is  introduced  from  the 
random  number  generator  wliicli  produces  erroneous  results  in  tJie  .simulation.  Tlie 
SVS-Fortran  random  number  generator  met  the  test  criteria  and  proved  to  be  one  of  the 
bettor  generators.  However,  this  generator  has  one  significant  draw  back.  'I'lie  random 
sequences  begin  at  one  of  two  dilVerent  values  depending  on  whether  the  seed  value  is 
positi\c  or  negatiaw.  'fhus  it  is  critical  llial  tlic  random  numbers  be  called  continuously 
in  a  loop  to  avoid  restarting  the  sequence,  thereby  introducing  unwanted  correlation. 

2.  Fast  Fourier  Transform 

The  most  repeatedly  used  algorithm  throughout  the  numerical  simulation  was 
the  fast  fouricr  transform.  Therefore,  it  was  necessary  to  use  the  most  elTicicnt  algorithm 
available.  The  Seitceh  Scientific  math  package  provided  several  options,  with  subroi.'inc 
II'TJC  .  best  suited  for  this  numerical  simulation.  This  subroutine  uses  a  complex  array 
input.  The  otltei  FFT  considered  was  a  routine  coded  by  Dr.  W'alters  which  he  rccei\  ed 
in  a  demonstration  package  provided  by  infotek.  This  TFT  utilizes  real  arrays  and  will 
henceforth  be  refered  to  as  subroutine  FFT . 

Faeh  subroutine  was  timed  for  various  configurations.  Specific  timing  results 
arc  contained  in  .Appendix  B.  Some  general  results,  however,  are  that  subioutinc 
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FFT2C  was  faster  for  one-dimensional  cases,  with  subroutine  FFT  faster  for  the  two- 
dimensional  cases.  The  decline  in  efficiency  of  FFT2C  can  be  attributed  to  the  extra 
coding  required  to  convert  between  real  arrays  and  complex  arrays.  Subroutine  FFF 
was  selected  over  FFT2C  since  the  simulation  utilized  the  two-dimensional  form  of  the 
FFT. 

Other  techniques  which  were  employed  to  increase  the  eHiciency  included  in¬ 
stalling  a  Weitek  coprocessor  in  the  Compaq.  This  reduced  the  processing  time  to  ap¬ 
proximately  one-third  that  of  the  original  time.  The  use  of  common  blocks  vice 
dimension  statements  further  decreased  processing  lime  by  5°  o.  Finally,  a  portion  of  the 
FFT  algorithm  was  modified  from  wim  =  sin(u/;g)  ,  to  wim  =  ^  1  —  wrc-  ,  with  a  negli¬ 
gible  decline  in  efficiency  by  0.0  F‘o. 

Because  of  the  discrete  nature  of  the  simulation,  there  exists  problems  and  lim¬ 
itations  associated  with  implementating  the  FFT.  One  such  problem  was  that  of  clas¬ 
sical  edge  diffraction  associated  with  the  phase  or  amplitude  discontinuities  at  the  edges 
of  the  finite  screen.  .As  the  propagation  distance  z  increases,  the  edge  diffraction  spreads 
toward  the  center  of  the  screen  making  more  and  more  of  the  diffraction  pattern  erro¬ 
neous.  Buckley  [Ref  1 1]  defines  the  distance  from  the  ends  of  the  screens  where  the  edge 
dilfraction  is  important  as 

D(rl  =  2^  ~  +  (>-^1 

where  z  is  the  piopagation  distance  and  </>„  is  the  root  mean  square  phase  deviation  im¬ 
posed  on  the  wa\e  by  the  screen.  The  severity  of  edge  effects,  however,  is  reduced  by 
the  aliasing  introduced  in  the  FFT  implementation.  Aliasing  transforms  the  linear 
screen  to  a  'circular'  one  with  the  last  point  associated  with  the  first.  This  results  in  a 
continuity  of  phase  and  amplitude  at  the  edges  of  the  screen  [Ref  1  !j. 

The  most  important  limitation  was  the  finite  spatial  range  imposed  by  the 
maximum  available  grid  size.  Hiis  places  a  constraint  on  the  asailable  range  of  fre¬ 
quencies  used  in  the  FFT  from  the  lowest  given  by  /„„„=1/Z.  ,  to  the  highest, 
=  njZL.  where  L  is  the  grid  length  and  n  is  the  number  of  grid  points.  Figure  ?  on 
page  20  illustrates  this  setup.  The  FF  l'  provides  a  least  squares  fit  of  sine  and  cosine 
fuctions  to  the  perturbed  wavefront  phases.  However,  this  method  prevents  an  accurate 
representation  of  the  wavefront  at  low  frequencies  for  a  Kolmogorov  x  "  ’’  power  spec¬ 
trum.  To  alleviate  this  problem,  a  subaperture  was  superimposed  at  the  center  of  the 
grid,  file  subaperture  helps  low  frequencies,  wiiicli  contain  a  large  portion  of  the  am- 
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plitude,  but  hurts  the  high  frequencies  by  limiting  the  inner  scale.  There  exists  some  high 
frequency  edge  effects,  hon'ever,  these  are  minimal. 

C.  PROCEDURE 

Figure  2  provides  a  summarx’  of  the  coding  contained  in  Appendix  C.  This  con¬ 
ceptual  diagram  illustrates  an  overview  of  the  procedural  steps  of  the  Fraunhofer  prop¬ 
agation  algorithm. 

1.  Input  Parameters 

Since  this  step  is  straightforward,  extensive  discussion  is  not  required.  However, 
it  is  important  to  note  that  the  input  parameters  were  both  fixed  and  variable.  The  array 
size  and  filter  value  were  fixed  quantities,  but  the  subaperture  size,  seed  value,  Q  value, 
and  propagation  distance,  took  on  dificrent  values.  The  actual  variable  names  are  doc¬ 
umented  at  the  beginning  of  the  code. 

2.  Aperture  Mutual  Coherence  Function 

The  second  step  in  this  procedure  calculates  the  aperture  MCF.  Subroutine 
MCF  does  this.  In  this  subroutine  the  initial  wavefront  was  represented  in  the  computer 
as  an  L  x  L  square  array  of  complex  numbers.  Centered  within  this  complex  electric 
lield  was  a  subapcrturc.  The  initial  complex  electric  field  had  a  value  of  zero,  every¬ 
where.  except  for  the  real  part  of  the  subaperturc,  which  had  the  value  of  one. 
Figure  3  illustrates  this.  With  the  electric  field  created,  it  was  direct  Fourier  transformed 
(DFT)  by  subroutine  DITIFT  .  The  intensity  of  the  electric  field  was  calculated,  and 
then  iir.crse  Fourier  transformed  (IFT).  yielding  the  aperture  .MCF. 

3.  Planar  Electric  Field 

The  complex  electric  field  was  created  by  the  same  method  prescubed  in  sub¬ 
routine  MCF  .  It  is  important  to  realize  that  the  concept  of  aperture  size  was  used  in 
two  ways.  One  way  corresponds  to  the  simulated  aperture  while  the  other  pertains  to 
an  aperture  with  physical  dimensions.  The  simulated  aperture  is  actually  a  matrix  or 
grid  which  directly  reflects  the  dimensioning  size.  For  example,  the  simulated  L  x  L 
complex  electric  field  was  actually  a  two-dimensional  matrix  corresponding  to  a  256  x 
256  two-dimensional  array.  From  the  input  parameters,  the  simulated  subaperture  takes 
on  grid  sizes  ranging  from  an  8  x  8  to  200  x  2U0  square  matrix. 

The  second  referencing  to  an  aperture  refers  to  an  aperture  with  actual  physical 
units.  The  physical  subaperture  was  assigned  a  value  of  0.3125  meters,  which  remains 
fixed  regardless  of  the  simulated  subaperturc  size.  The  length  L  of  the  complex  electric 
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Figure  2.  Conceptual  Diagram  of  the  Simulation  Coding. 
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field,  on  the  other  hand,  had  variable  physical  lengths  depending  on  the  simulated 
supaperture  size.  The  physical  value  of  L  was  calculated  from 


L-ntx^-. —  [meters].  (5.>) 

tsize 

where  m  is  the  physical  subaperture  length  in  meters,  nr  is  the  integer  value  of  the  array 
dimensioning,  and  isize  is  the  integer  value  for  the  simulated  subaperture. 

4.  Phase  Screens 

a.  Generation 

Phase  screens  were  created  in  subroutine  GGAUS  by  constructing  a  L  x  L 
matrix  of  Gaussian  distributed  random  numbers.  Each  position  was  assigned  an  inde¬ 
pendent  random  number  n,  thereby,  requiring  n^  random  numbers.  Since  the 
SVS-Fortran  random  number  generator  was  uniformly  distributed,  an  algorithm  pro¬ 
vided  by  Knuth  (Ref.  10)  tramsformed  the  distribution  into  a  Gaussian  one.  Two  inde¬ 
pendent.  two-dimensional  real  arrays  called  phaser  and  phase:  were  created,  which 
represent  the  real  and  imaginaiy-  components  of  a  two-dimensional  complex  phase 
screen.  In  this  algorithm,  the  imaginary  part  of  the  phase  screen  was  set  to  zero. 

The  domain  in  which  the  phase  screens  are  created  is  arbitraiy.  however, 
filtering  was  done  in  the  Fourier  plane.  Creating  the  phase  screen  in  frequency  space 
vice  real  space,  reduces  the  requirement  for  an  additional  FFT  when  transforming  from 
real  space  to  frequency  space.  Martin  and  Flatte  [Ref  7]  proposed  this  method,  but  it 
creates  dilficulties  in  absolute  normalization.  Further  discussion  is  presented  in  Chapter 
four.  This  simulation,  on  the  other  hand,  generated  the  random  Gaussian  phase  screen 
in  real  space.  This  in  turn  was  DFT'd  to  frequency  space  where  the  complex  phase 
screen  was  filtered  and  then  the  filtered  phase  screen  was  IFT'd. 

b.  Filtering 

Phase  screens  were  filtered  to  obtain  the  correct  power  law  form.  The  fil¬ 
tering  function  used  in  subroutine  FLTR  was, 

=  2nk^dj^O„.  (56) 

where  S,  is  the  slab  thichness  and  =  0.0.13Qk'" ,  which  gives  the  relationship  be¬ 
tween  the  phase  spectrum  and  refractive-index  spectrum  [Ref  2.  pp.  101-102.].  Filtering 
was  accomplished  by  multiplying  each  phase  screen  spectrum  by  the  square  root  of 
equation  (56).  The  correct  filtering  method  requires  circular  frequency  filtering  instead 
of  a  linear  one.  because  of  two-dimensional  isotropy.  The  correct  form  is 
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K  =  ^  K-  +  k;  ,  which  is  radial  ever\Avhcre  except  at  the  origin,  where  it  is  zero  [Ref.  12j. 
This  was  reflected  in  the  simulation  by  setting  the  position  (1.1)  equal  to  zero  in  each 
two-dimensional  real  array  which  makes  up  the  complex  phase  screen. 

It  is  important  to  realize  that  although  the  filtering  concept  is  simple  and 
straightforward,  the  actual  implementation  is  not.  The  difficulty  arises  from  the  sym¬ 
metry  properties  of  the  digitally  filtered  phase  screen. 

A  one-dimensional  case  offers  a  simple  illustration  of  this  concept.  The 
FFT  of  a  complex  function  which  contains  only  real  components,  results  in  a  symmetric 
function  in  frequency  space  about  the  Nyquist  frequency  for  the  real  components,  and 
an  anti-symmetric  function  for  the  imaginaries.  With  these  symmetries  present  in  the 
frequency  domain,  the  correct  implementation  of  the  filtering  is  to  mimick  these  sym¬ 
metries.  Thus  a  "folding"  technique  about  the  Nyquist  frequency  was  required.  How¬ 
ever.  when  this  concept  was  extended  to  a  two-dimensional  case,  as  in  the  phase  screen 
in  the  simulation,  the  symmetries  imposed  by  the  FFT,  were  no  longer  apparent  in  the 
frequency  domain.  To  simplify  the  symmetn."  requirements,  a  real  phase  screen  was 
used.  The  real  and  imaginary  spectral  components  were  filtered  by  folding  about  the 
Nyquist  frequency.  It  was  suggested  by  both  Brigham  [Ref  13]  and  Martin  and  Flatte 
that  a  complex  phase  screen  containing  both  real  and  imaginary  components,  yields  two 
entirely  distitict  phase  screens.  However,  nothing  was  procided  to  support  this  hypoth¬ 
esis. 

.Anotlior  important  consequence  of  filtering  resides  in  the  units.  The  phase 
screens  were  filtered  in  k  units,  but  the  FFT  algoritiim  operates  in  frequency  units. 
Therelbre.  it  was  necessary  to  make  a  change  of  variables  prior  to  applying  the  ITT. 
The  relationship  used  for  the  change  of  variables  is  k  =  2nv  . 

The  paper  by  Martin  and  Flatte  [Ref  7].  specifies  an  additional  nermaliza- 
tion  factor  of  A:'  in  equation  (56).  where  A,  =  where  N  is  the  number  of  grid 

points  and  A  is  the  sampling  interval.  Martin  and  Flatte  provided  no  explanation  for 
the  additional  A:'  in  the  filtering.  It  was  not  included  in  the  filtering  code. 

5.  Impleinenlation 

After  the  filtered  phase  screen  was  IFT'd  into  real  space,  it  was  introduced  into 
the  code  as  a  phase  screen  and  multiplied  with  the  electric  field.  1  he  array  phaser,  which 
contains  tiie  desired  phase  field,  takes  the  form  of  the  Rytov  approximation,  .  in  the 
extended  1  luygens-Frescnel  integral.  This  form  assumes  that  only  phase  perturbations 
and  not  amplitude  variations,  occur. 


6.  Propagation  Methods 

As  indicated  in  previous  sections,  the  Huygens-Fresnel  technique  was  used  to 
simulate  the  propagation  of  light.  This  was  accomplished  by  applying  the  TFT  to  the 
perturbed  electric  field.  Both  the  "far-field"  and  "near-field"  propagation  methods  were 
considered. 

a.  Fraunhofer 

Of  the  two  different  propagation  methods,  the  single  screen  Fraunhofer 
propagation  is  by  far  the  simplest  technique  to  implement  and  the  one  implemented  in 
this  thesis.  The  uniform  coherent  plane  wave  at  the  aperture  was  FFT'd  yielding  the 
desired  diffraction  pattern.  Looking  more  closely,  one  can  see  that  under  certain  cir¬ 
cumstances,  Fraunhofer  propagation  is  just  a  special  case  of  the  long  distance  convo¬ 
lution  form  of  Fresnel  propagation  given  by  equation  (46).  There  exists  two  situations 
when  this  occurs.  One  is  a  "far-field"  case  for  large  distances,  where  the  point  of  obser¬ 
vation  is  at  infinity.  1  he  other  case,  is  when  a  spherical  curvature  is  placed  on  the  wave 
at  the  aperture.  This  curvature  cancels  the  quadratic  phase  factor  in  the  Fresnel  form, 
at  the  focal  point. 

Both  Fresnel  and  Fraunhofer  algorithms  were  needed  for  propagation.  It 
was  essential  to  verify  and  validate  each  case  before  building  on  the  pre-existing  codes. 
The  Fraunhofer  algorithm  provided  the  basis  of  this  simulation.  Since  the  Fraunhofer 
dilT'raction  pattern  is  well-known,  it  provided  a  means  to  verify  the  existing  code  by 
comparing  the  simulated  dilTraction  pattern  with  theoretical  results.  Figure  4  and  Fig¬ 
ure  5  illustrate  one  three-dimensional  quadrant  of  the  dilTraction  pattern  of  an  unper¬ 
turbed  electric  field,  for  two  diflcrcnt  subaperture  sizes.  \\'hile  Figure  6  illustrates  one 
thrce-din:ensional  quadrant  ol  a  perturbed  electric  field  dilTraction  pattern  for  a  16  \  16 
subaperture. 

b.  Fresnel 

Although  the  Fresnel  propagation  forms  were  implemented  but  not  tested 
in  this  thesis,  discussion  is  warranted  since  multi-screen  Fresnel  propagation  is  predom¬ 
inantly  used  in  thermal  blooming  and  multiple  scattering  scenarios.  Propagation 
through  the  turbulent  boundary  layer  also  requires  Fresnel  propagation  codes.  Two 
diflerent  forms  are  used  for  Fresnel  propagation,  which  apply  a  straightforward  FFT  to 
evaluate  them.  These  two  forms  however,  do  not  allow  for  a  variable  receiving  array 
size.  In  addition:  choosing  the  correct  number  of  sample  points  is  essential.  An  elTective 
approach  used  to  resolve  these  problems  is  the  Fresnel  number. 
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Figure  5. 


Three  Dimensional  Diffraction  Pattern  of  a  16  x  16  Subaperture. 


Figure  6, 


Perturbed  Diffraction  Pattern  of  a  16  x  16  Subaperture. 


The  Fresnel  number  is 


where  A  r  is  the  receiving  aperture  size,  A  p  is  the  transmitting  aperture  size.  ).  is  the 
wavelength  of  light,  and  z  is  the  propagation  distance.  Implementation  of  the  correct 
Fresnel  form  depends  on  the  Fresnel  number.  When  the  Fresnel  number  is  smaller  than 
the  number  of  grid  points  in  the  field  length,  that  is,  ,  the  long  distance  propa¬ 

gation  algorithm  is  used.  Conversely,  when  A'>A'}  ,  the  short  distance  code  is 
implemented. [Ref  14  ] 

The  long  distance  propagation  code  uses  the  convolution  form  of  the 
diffraction  integral.  Implementation  requires  placing  a  curvature  on  the  electric  field 
wavefront  at  the  aperture.  The  subroutine  called  quadl  ,  docs  this.  Multiplying  the 
phase  screen  with  the  electric  field,  gives  a  perturbed  electric  field  that  propagates  the 
entire  distance  by  means  of  one  FFT.  In  the  Fourier  plane,  the  quadratic  phase  factor 
called  quadl  scales  the  electric  field. 

The  transfer  function  form,  on  the  other  hand,  is  suited  for  short  distances. 
In  this  case,  the  entire  propagation  distance  is  divided  into  equally  spaced  slabs.  Two 
FFTs  arc  required  to  propagate  the  distance  of  each  slab.  This  is  accomplished  by  the 
following  method: 

1.  Mesh  the  electric  field  and  phase  screen. 

2.  .Apply  the  direct  FFT. 

3.  Multiply  the  field  by  the  propagation  transfer  function  subroutine  called  irmfr  . 

4.  Apply  the  invei'-e  FFT. 

1  his  procedure  is  repeated  for  each  phase  screen  until  the  observing  plane  is  reached. 
[Ref  7) 

Both  forms  of  Fresnel  propagation  were  included  in  the  simulation  code 
contained  in  Appcndi.x  C,  however,  neither  form  of  Fresnel  propagation  was  e.vercised. 

7.  Atmospheric  Mutual  Coherence  Function 

The  final  step  in  this  simulation  calculated  the  atmospheric  .MCF.  The  same 
procedure  used  to  calculate  the  aperture  MCF,  was  applied  to  the  perturbed  electric 
field,  with  one  exception.  The  diflerence  is  that  division  of  the  composite  atmospheric- 
optical  MCF  by  the  aperture  MCF.  gives  the  atmospheric  .MCF. 


IV.  RESULTS 


The  main  objectives  of  this  thesis  were  to  demonstrate  that  the  simulation  gives  re¬ 
sults  that  provide  insight  to  weak  scattering  wave  propagation  and  examine  the  accuracy 
of  the  simulation  code.  The  limitations  imposed  by  the  computer  mechanics  and  the  ' 
actual  method  of  implementing  turbulence  theory  are  discussed.  Specific  areas  of  inter¬ 
est  include  the  filter  function,  MCF,  and  saturation  efi'ects. 

A.  FILTERING 

As  previously  mentioned  in  chapter  three,  Martin  and  Flatte  proposed  a  vers'  dif¬ 
ferent  approach  to  creating  the  random  Gaussian  phase  screen.  Their  method  suggests 
creating  the  phase  screen  in  frequency  space,  vice  real  space,  to  reduce  the  required 
number  of  FFT's  from  two  to  one.  Since  the  domain  in  which  the  phase  screen  is  gen¬ 
erated.  is  arbitrary,  this  approach  seems  plausible.  However,  this  method  proved  to  he 
awkward.  As  the  subaperture  size  was  successively  doubled,  it  was  necessar}’  to  increase 
the  strength  of  turbulence,  Q,  by  a  factor  of  ten  each  time,  in  order  to  produce  the 
identical  level  of  turbulence  as  in  the  previous  phase  screen.  Additionally,  since  the 
phase  screen  was  created  in  frequency  space,  and  no  symmetries  were  present,  the  fil¬ 
tering  teclmique  did  not  reflect  any  folding  about  the  Nyquist  frequency.  It  is  not  clear 
that  aliasing  was  accounted  for  in  the  NFartin  and  Flatte  algorithm.  Hence,  the  IFT  of 
the  phase  screen  appears  to  result  in  a  statistically  incorrect  phase  screen.  Additionally, 
it  is  not  clear  hon-.  with  one  FFT.  Martin  and  Flatte  handled  the  2r.  and  1  N  normril- 
ization  reqiiircmeitts.  N\'ith  a  round  trip  of  FFT's,  the  normalization  problems  are  au¬ 
tomatically  handled.  The  weakly  filtered  phase  screen,  in  turn,  led  to  MCF  curves  which 
gros''Iy  overestiinatcd  the  coherence  length.  /Aj-  These  problems  obtained  from  .Martin 
and  Flatte  '  s  approach  to  simulating  turbulence  led  to  the  current  coding  which  created 
one  phase  screen  in  real  space  and  applied  a  folding  about  the  Nyquist  frequency,  in  the 
filtering,  to  account  for  the  symmetries  introduced  by  the  FFT's. 

Kolmogorov  theory  of  atmospheric  turbulence  predicts  that  the  graph  of  the  phase 
screen  power  spectral  density  versus  k  yield  a  slope  of —11/3.  Figure  7  shows  that 
equation  (56)  produced  a  filtered  phase  screen  that  reflects  the  —11/3  slope,  as  well  as, 
the  correct  folding  technique.  Identical  slopes  were  expected  for  all  subapertures,  as  well 
as.  all  possible  angles  which  reflect  the  circular  filtering.  Other  subaperture  profiles  show 
a  consistent  slope  value  of -11/3.  Figure  S  corresponds  to  a  32  \  32  subaperture  at  a 


45  degree  angle  while  Figure  9  is  representative  of  a  64  x  64  subapcriure  at  a  90  degree 
angle.  Isotropy  is  apparent  in  that  the  circular  filtering  was  implemented  correctly 
within  the  euor  introdui-cd  by  the  randomness  in  the  screen  and  the  ability  to  linearly 
fit  a  line  through  the  data  points. 

B.  MCF 

The  MCF  of  the  electric  field  provides  one  method  of  analyzing  the  accuracy  of  the 
simulation.  To  verify  that  the  MCF  was  correctly  computed,  the  aperture  .MCF  of  an 
unperturbed  electric  field  was  calculated.  This  was  easi]>-  accomplished  since  the  image 
intensity  and  aperture  .MCF  are  transform  pairs  and  are  analytical  for  simple  square  and 
circular  apertures.  The  aperture  .MCF  is  just  the  autocorrelation  of  the  aperture  func¬ 
tion  which  is  evaluated  by  calculating  the  area  of  overlap  of  two  identical  apertures  as 
they  are  moved  laterally  apart.  For  a  square  subaperture,  the  autocorrelation  yields  the 
triangle  function,  with  maximum  value  of  1.0  and  ntinimum  value  of  0.0  corresponding 
directly  to  the  subaperture  size.  Figure  10  illustrates  the  MCF.  or  autocorrelation,  of 
an  8  X  S  and  16  x  16  square  subaperture. 

The  .MCF  corresponding  to  the  atmospheric  turbulence  was  deternrined  by  dividing 
the  .MCF  of  the  perturbed  electric  field  by  the  aperture  .MCF.  The  value  of  the  spa¬ 
tial  coherence  length,  was  determined  from  the  turbulence  MCF  curve.  Figure  11  il¬ 
lustrates  a  simulated  value  of  .'.20  nun  for  a  64  x  64  subaperture  with  C;  = 

The  theoretical  value  calculated  from  equation  (40)  yields  p,  =  2.41/n/n.  Although  the 
64  X  64  subaperture  gave  accurate  results  other  subapcriure  configurations  did  not. 
When  0:  increased,  the  .MCF  curves  fell  olT  rapidly  towards  zero  for  all  subapertures. 
Figure  12  illustrates  this  for  a  64  x  64  subaperture.  .Ml  subaperture  configurations  were 
run  for  two  dilTerenct  Q.  values.  The  simulated  po  values  were  plotted  against  the  cor¬ 
responding  subaperturcs  for  each  case.  Figure  13  reflects  the  Q  =  l.xTO  'Vo  values,  and 
Figure  14  is  for  C,  =  l.vlO'''.  The  desired  trend  is  for  the  simulated  Po  values  to  approach 
the  theoretical  value  as  the  subaperture  size  increased.  This  trend  is  visible  in 
Figure  13  where  Q  =  l.vlO''^  and  p,  =  2.4mm.  Ifowevcr,  Figure  14  shows  continuously 
decreasing  p^  values  past  the  theoretical  value  of  9.6mm.  The  results  of  figures  thirteen 
and  fourteen  point  to  a  problem  that  may  involve  edge  effects  or  aliasing,  resulting  from 
undersampling.  The  p,j  values,  which  were  on  the  order  of  millimeters,  were  smaller  than 
the  tens  of  centimeter  distances  of  the  subaperture  mesh  size. 

In  an  attempt  to  pinpoint  the  problem,  the  simulation  code  was  changed  to  increase 
the  number  of  frequencies  in  the  subaperture  by  using  a  512  x  512  array.  The  results 
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Figure  7.  Filtered  Phase  Screen  of  a  16  x  16  Subaperture. 
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Filtered  Phase  Screen  of  a  32  ,x  32  Subaperture. 
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Figure  1 1.  Mutual  Coherence  Function  of  a  64  x  64  Subaperture. 
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were  identical  to  those  for  the  256  x  256  dimensioning,  within  the  arithmetic  error  of  the 
algorithm.  The  information  provided  by  the  512  dimensioning  was  that  the  inaccuracy 
of  the  code  was  not  due  to  an  inner  scale  problem,  however,  a  low  frequency,  outer  scale 
problem  may  still  exist. 

C.  SATURATION 

The  last  area  of  investigation  was  whether  or  not  the  simulation  predicts  the  satu¬ 
ration  of  intensity  for  an  extended  medium  modeled  by  a  single  phase  screen  realization. 
Saturation  is  generally  considered  to  be  caused  by  multiple  scattering,  a  scattered  wave 
interfering  with  a  distorted  wave.  One  would  expect  from  the  Rytov  approximation  of 
turbulence  theory,  that  saturation  will  occur  even  in  Fraunhofer  propagation.  This  as¬ 
sumption  stems  from  the  representation  of  turbulence  in  the  form,  e'®,  which  has  a 
magnitude  bounded  between  plus  and  minus  one.  Figure  15  illustrates  that  the  nor¬ 
malized  intensity  variance  saturates  with  increasing  turbulence.  Theoretically,  a  nor¬ 
malized  variance  of  one  is  expected  for  Rayleigh  statistics.  However,  it  is  premature  to 
assume  that  saturation  is  inherent  in  Fraunhofer  cases,  until  the  MCF  results  arc  veri¬ 
fied. 
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Figure  15.  Intensity  Variance  Saturation  of  a  16  x  16  Subaperture. 
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V.  CONCLUSIONS  AND  RECOMMENDATIONS 

This  thesis  simulated  the  propagation  of  plane  waves  through  an  extended  two- 
dimensional  random  media  using  a  single  phase  screen  technique.  The  atmospheric 
turbulence  had  a  Kolmogorov  pure  power  law,  while  the  propagation  was  strictly 
Fraunhofer.  Limitations  of  its  applicability  were,  primarily,  the  finite  spatial  range  im¬ 
posed  by  the  available  grid  size.  Diffraction  patterns,  correlation  functions  and  intensity- 
variance  saturation  at  the  observation  plane,  were  investigated. 

The  results  provided  by  the  simulation  suggest  general  agreement  with  the  turbu¬ 
lence  theory.  Saturation  for  weak  scattering  was  supported  by  this  model.  The  .\1CF 
curves,  although  not  completely  correct,  provided  insight  to  the  theory  and  illustrated 
problems  which  are  still  present  in  the  current  coding.  Some  specific  problems  include, 
potential  errors  in  the  implementation  of  the  filtering  technique  from  edge  effects,  alias¬ 
ing  and  inner  scale  problems,  or  from  incorrect  normalization. 

.Any  further  research  on  this  topic  should  begin  with  resolving  the  inaccuracies  still 
present  in  the  simulation  coding.  Several  possible  reasons  were  presented,  however,  an 
error  in  the  filtering  of  the  phase  screen  seems  to  be  the  most  likely  cause.  Further 
testing  can  be  conducted  on  the  phase  screen,  to  include  calculating  the  phase  screen 
variance,  as  well  as.  its  structure  function  Z)„.  By  comparing  the  simulated  phase  screen 
with  the  theoretical  structure  function  for  D^.  this  technique  will  verify  whether  or  not 
the  simulated  phase  screen  accurately  represents  turbulence. 

An  assumption  was  made,  that  aliasing  was  not  a  problem,  since  the  problems  as¬ 
sociated  with  undersampling  were  not  apparent.  -Aliasing  can  he  tested  by  using  liner 
grid  sizes  and  observing  changes  induced  by  the  higher  spatial  frequencies. 

-After  the  coding  is  working  correctly,  both  the  convolution  and  transfer  function 
form  of  Fresnel  propagation  can  be  implemented  and  exercised.  In  addition,  the  array 
sizes  should  be  modified  to  incorporate  a  dimension  of  512  or  1024.  Finally,  during  the 
phase  screen  generation,  phasei  should  be  filled  with  random  numbers  to  give  two  usable 
phase  screens.  It  is  not  obvious  whether  two  independent  phase  screens  will  be 
produced,  or  whether  the  total  energy  will  be  distributed  among  the  two  phase  screens. 
Therefore,  testing  should  be  conducted  to  ensure  that  each  usable  phase  screen  possesses 
the  correct  statistics. 
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APPENDIX  A.  RANDOM  NUMBER  GENERATOR  TEST  DATA 


The  following  table  provides  the  statistical  results  of  the  five  empirical  tests  run  on 
the  random  number  generator,  Ran^Ij.  The  values  were  determined  from  the  Chi- 
Square  table  in  Bevington  [Ref  15  ]. 

The  results  of  the  Lagged  Product  test  correspond  respectively  to  the  theoretical 
mean,  (jlj  ,  calculated  mean,  fx  ,  theoretical  standard  deviation,  a-r  ,  and  calculated 
standard  deviation,  a.  A  lag  of  three  was  tested. 


Table  I.  RANDOM  NUMBER  GENERATOR  TEST  DATA 


TEST 

DEGREES  OF 
FREEDOM 

PROBABILITY 

Frequency 

9 

4.4 

SS.3% 

Serial 

20 

17.1 

65.0"  0 

l.’p  and  Down 

8 

4.1 

84.5% 

Above  and  Below  the 
Mean 

12 

Bi 

90.5% 

TEST 

l^r 

a 

Lagged  Product 

0. 2.5(1 

0.259 

0.1 00 

0.U92 
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APPENDIX  B.  FFT  TIME  SERIES  DATA 


This  appendix  contains  the  results  of  the  comparison  between  the  two  FFT  sub¬ 
routines,  FFT2C  and  subroutine  FFT  .  A  function  routine  called  SCXDS,  was  imple¬ 
mented  in  the  code  to  provide  accurate  time  measurements. 


Table  2.  FFT  TIME  SERIES  DATA 


COMPUTER  SETUP 

FFT 

FFT2C 

1  DIMENSION 

TIME(sec) 

TlME(sec) 

2>5  with  20-MH2  80387 

2.64 

2.64 

2‘'^  with  20-MHz  S03S7 

25.  .^4 

25.26 

2'-  with  20-MHz  S03S7 

54.05 

53.44 

2  DIMENSION 

Tl,ME(sec) 

TlMEisec) 

12S  X  12S  with  20-MHz  S03S7 

6. 86 

7.08 

256  X  256  with  20-MHz  S0.vS7 

31.30 

}2M 

256  \  256  with  W'eitek 

12.0S 

12.14 

256  X  256  with  Weitok  and  common 
block 

11.48 

11.56 

512  X  512  with  \\'citek  and  comrnon 
block 

40.54 

52.40 
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APPENDIX  C.  SIMULATION  CODE 


The  simulation  code  contained  in  this  appendix  incorporates  Fraunhofer  and  both 
forms  of  Fresnel  propagation.  This  thesis  only  exercised  the  Fraunhofer  propagation. 


THIS  CODE  PROVIDES  A  QUALITATIVE  VIEW  OF  BOTH  FRAUNHOFER  AND 
FRESNEL  DIFFRACTION  BY  OBSERVING  THE  PERTURBATION  IMPOSED  ON  A 
MONOCHROMATIC  PLANE  WAVE  PROPAGATING  THROUGH  A  TURDULENT 
MEDIUM.  THE  TURBULENT  MEDIUM  IS  INTRODUCED  IN  THE  FORM  OF  A 
STOCHASTIC  PHASE  SCREEN.  PROPAGATION  OF  THE  ELECTRIC  FIELD  IS 
ACCOMPLISHED  THROUGH  FFT'S. 


GLOSSARY  OF  VARIABLE  NAMES: 

1.  ARRAYS; 

RE  -  ONE  DIMENSION  REAL  ARRAY  OF  DIMENSION  NR,  WHICH  IS  USED  TO 
MANIPULATE  THE  REAL  PART  OF  THE  COMPLEX  ELECTRIC  FIELD  IN 
THE  FFT  SUBROUTINE.  THIS  ARRAY  IS  REPEATEDLY  USED  THROUGHOUT 
THE  CODE. 

RIM-  ONE  DIMENSION  REAL  ARRAY  OF  DIMENSION  NR,  WHICH  IS  USED 
TO  MANIPULATE  THE  IMAGINARY  PART  OF  THE  COMPLEX  ELECTRIC 
FIELD  IN  THE  FFT  SUBROUTINE.  THIS  ARRAY  IS  REPEATEDLY  USED 
THROUGHOUT  THE  CODE. 

FIELDR  -  TtCO  DIMENSION  REAL  ARRAY  OF  DIMENSION  NR  X  NR 

CONTAINING  THE  REAL  PART  OF  THE  COMPLEX  ELECTRIC 

FIELD.  THIS  ARRAY  IS  REPEATEDLY  USED  THROUGHOUT  THE  CODE 

FIELDI  -  TVv^O  DIMENSION  REAL  ARRAY  OF  DIMENSION  NR  X  NR 

CONTAINING  THE  IMAGINARY  PART  OF  THE  COMPLEX  ELECTRIC 
FIELD.  THIS  ARRAY  IS  REPEATEDLY  USED  THROUGHOUT  THE  CODE 

FILL  -  TVi'O  DIMENSION  REAL  ARRAY  OF  DIMENSION  NR  X  NR  USED  AS  A 

DUMMY  ARRAY  IN  THE  IFT  OF  THE  POWER  SPECTRUM  WHICH  YIELDS 
THE  MCF. 

FIELDM  -  TW'O  DIMENSION  ARRAY  OF  DIMENSION  NR  X  NR  REPRESENTING 
THE  MAGNITUDE  OF  THE  PERTURBED  ELECTRIC  FIELD. 

FMCF  -  TV*’0  DIMENSION  ARRAY  OF  DIMENSION  NR  X  NR  REPRESENTING  THE 
INTENSITY  OF  THE  PERTURBED  ELECTRIC  FIELD.  THIS  ARRAY  IS 
USED  IN  DETERMINING  THE  MCF. 

FNORM  -  TW'O  DIMENSION  ARRAY  OF  DIMENSION  NR  X  NR  REPRESENTING  THE 
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INTENSITY  OF  THE  UNPERTURBED  ELECTRIC  FIELD.  THIS  ARRAY 
IS  ALSO  USED  IN  DETERMINING  THE  MCF. 

PHASER  -  TVv’O  DIMENSION  REAL  ARRAY  OF  DIMENSION  NR  X  NR  CONTAINING 
THE  REAL  PART  OF  THE  RANDOM  COMPLEX  PHASE  SCREEN. 

PHASEI  -  TWO  DIMENSION  REAL  ARRAY  OF  DIMENSION  NR  X  NR  CONTAINING 
THE  IMAGINARY  PART  OF  THE  RANDOM  COMPLEX  PHASE  SCREEN. 

FMAG  •  ONE  DIMENSION  SLICE  OF  THE  PERTURBED  MCF  USED  IN  THE 
GRAPHICS  ROUTINE.  REAL  ARRAY. 

DIST  -  ONE  DIMENSION  REAL  ARRAY  REPRESENTING  THE  PIXELS 
CORRESPONDING  TO  THE  VALUES  IN  THE  ARRAY  FMAG. 

VARAIBLES: 

NR  -  DIMENSION  OF  THE  ARRAYS  EXACTLY  AS  SPECIFIED  IN  THE  DIMENSION 
STATE.MENTS  IN  THE  CALLING  PROGRAM.  INPUT  INTEGER.  INDICE. 

N2  -  ONE  HALF  OF  NR.  INTEGER. 

M  -  POWER  OF  2  IN  THE  DIMENSIONING.  USED  IN  THE  FFT  SUBROUTINE. 

ISIZE  -  INTEGER  VALUE  CORRESPONDING  TO  A  PARTICULAR  CHOICE  FOR  AN 
APERTURE  SIZE.  INPUT  VARIABLE. 

NSIZE  -  INTEGER  VALUE  CORRESPONDING  TO  THE  ACTUAL  APERTURE  SIZE. 

DELMSH  -  REAL  VALUE  REPRESENTING  THE  SAMPLING  INTERVAL  FOR  A 
PARTICULAR  APERTURE  SIZE. 

SEED  -  REAL  INPUT  VARIABLE  USED  TO  BEGIN  A  RANDOM  SEQUE.NCE  OF 
NUMBERS  FOR  THE  SUBROUTINE  GGAUS. 

NYES  -  INPUT  INTEGER  VARIABLE  WHICH  SELECTS  WHETHER  TURBULENCE  IS 
INTRODUCED  IN  THE  CODE. 

FILTER  -  FIXED  INPUT  VALUE  USED  IN  THE  FILTERING  FUNCTION.  REAL. 

CN2  -  REAL  INPUT  VARIABLE  REPRESENTING  THE  INDEX-OF-REFRACTION 
STRUCTURE  PARAMETER.  DETERMINES  THE  AMOUNT  OF  TURBULENCE 
I.NTRODUCED  IN  THE  FILTERING  FUNCTION. 

DREC  -  REAL  INPUT  VARIABLE  REPRESENTING  THE  RECEIVING  FIELD  SIZE. 

DTRNS  -  FIXED  VALUE  FOR  THE  TRANSMITTING  FIELD  SIZE.  REAL. 

Z  -  REAL  INPUT  VARIABLE  REPRESENTING  THE  TOTAL  PROPAGATION 
DISTANCE  OF  THE  ELECTRIC  FIELD. 

DELX  -  REAL  VARIABLE  FOR  THE  PROPAGATION  DISTANCE  TO  EACH  SLAB. 

NUMSCR  -  INTEGER  INPUT  VARIABLE  USED  IN  NEAR-FIELD  FRESNEL 

PROPAGATION.  THIS  VARIABLE  CORRESPONDS  THE  THE  NUMBER 
OF  EQUALLY  SPACED  SLABS  WHICH  MAKE  UP  THE  TOTAL  DISTANCE. 
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c 

C  VVL  -  FIXED  VALUE  FOR  THE  WAVELENGTH. 

C 

C  ICNT  -  INTEGER  VALUE  WHICH  DETERMINES  WHETHER  OR  NOT  THE  FIELD  HAS 
C  PROPAGATED  THE  TOTAL  DISTANCE  Z.  USED  IN  NEAR-FIELD 

C  FRESNEL  PROPAGATION. 

C 

C  PI  -  VALUE  OF  PI. 

C 

C  TPI  -  TWICE  THE  VALUE  OF  PI. 

C 

C  MODE  -  REAL  VALUE  WHICH  DETERMINES  THE  FORM  OF  PROPAGATION. 

C 

C  SIGN  -  REAL  VALUE  EITHER  1.0  OR  -1.0  WHICH  DETERMINES  WHETHER  THE 
C  FFT  IS  DIRECT  OR  INDIRECT.  IT  ALSO  DETERMINES  WHETHER 

C  NORMALIZATION  OCCURS. 

C 

C  FLDM  -  MAXIMUM  VALUE  IN  THE  PERTURBED  MCF  ARRAY. 

C 

C  FMAX  -  MAXIMUM  VALUE  IN  THE  PERTURBED  ELECTRIC  FIELD  ARRAY. 

C 

C  GRAPHICS; 

C 

C  THE  FOLLOWING  VARIABLE  NAMES  ARE  EITHER  SPECIFIC  TO  THE  SVS 
C  GRAPHICS  ROUTINE  OR  USED  TO  MANIPULATE  DATA  FOR  GRAPHING: 

C 

C  NDEX,  IREG,  ANS,  GETC,  lUNITP.  lUNITV,  ISYMB,  ITNO,  MON,  NPRIN, 

C  MODE,  ONE,  TITLE,  IX,  lY,  XMAX,  ICOLOR,  INCR,  NTOT. 

C 

C 

C 

COMMON  /BLKl/  RE(256)  ,RI.M(256) 

COMMON  /BLK2/  FIELDR( 256 ,256) ,FIELDI( 256 , 256) 

COMMON  /BLK3/  PHASER( 256 ,256) ,PHASEI( 256 , 256) 

DIMENSION  NDEX(20) ,FIELDM(256,256) ,FMCF( 256 , 256) ,FNORM( 256 , 256) 
DIMENSION  FMAG( 130) ,DIST( 130) ,FILL(256,256) 

INTEGER*2  IREG(9) 

DOUBLE  PRECISION  PI 

DATA  NDEX/ 15, 15,7,7,8,8,14,14,5,5,4,4,2,2,3,3,1, 1,0,0/ 

DATA  RE/256-'-0.  0/ ,RIM/256*0.  0/ 

DATA  FMAG/130*0.  0/ ,DIST/ ISO^’-O.  0/ 

CHARACTER-'4  ONE 

CHARACTER*21  TITLE 

CHARACTER---2  ANS , GETC 

DATA  TITLE/ 'THE  SEED  VALUE  IS=  '/ 

DATA  ONE/ '  ' / 

DATA  IUNTTP/10/,IUNITV/20/ 

0PEN(4,FILE='  L.N.  DAT'  , STATUS= ' NEW '  ) 

ISYMB=22 

ITN0=5 

M0N=18 

PI=3. 141592653589792 

NPRIN=0 

MODE=0 
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999  CONTINUE 


THIS  SECTION  OF  THE  PROGRAM  SETS  UP  THE  INPUT  PARAMETERS  FOR  THE 
SIMULATION. 

WRITE(*,*)' 

ls'RITE( 'HELLO.  .  .  LET  US  BEGIN  THIS  SIMULATION  BY  E.NTERING  ’ 
WRITE(*,*)' SEVERAL  INPUT  PARAMETER  VALUES.' 

WRITE(*,*)' 

WRITEC*,*) 'THE  VARIABLE  WHICH  DIMENSIONS  THE  ARRAY  SIZE  IS  THE' 
WRIT£(*,*) 'FIRST  VALUE  TO  ENTER.  INPUT  THE  INTEGER  VALUE.  ' 

RE AD (*,*). NR 
WRITE (*,*)' 

WRITE(*,*) 'THE  SECOND  VARIABLE  OF  INTEREST  IS  "NSIZE".  THIS' 
WRITE(*,*)' VARIABLE  DIMENSIONS  THE  SIZE  OF  THE  PLANAR  ELECTRIC' 
WRITEC*,*) 'FIELD.  SELECT  ONE  OF  THE  FOLLOWING.  ' 

WRITE(*,*)'  1.  FOR  100  X  lOO’ 

WRITEC*,*)'  2.  FOR  64  X  64’ 

WRITEC*,*)'  3.  FOR  32  X  32' 

V.'RITEC*,*)'  4.  FOR  16  X  16’ 

WRITEC*,*)'  5.  FOR  8  X  8' 

WRITE r*,*)'  6.  FOR  4  X  4' 

WRITEC*,*)'  7.  FOR  2  X  2’ 

WRITEf*,*)' 

READC*,*'nSIZE 
IFCISIZE. EQ.  1)  THEN 
NSI2E=50 
DELMSH=.  0031 
INCR=3 

ELSEIFCISIZE.  EQ.  2)  THEN 
NSI2E=32 
DELMSH=.  C049 
I.NCR=1 

ELSEIFCISIZE.  EQ.  3)  THEN 
NSIZE=16 
DELMS1[=.  0098 
1NCR=2 

ELSEIFCISIZE.  EQ.  4)  THEN 
NSIZE=S 

d::lmsk=.  019531 

INCR=2 

ELSEIFCISIZE. EQ. 5)  THEN 
NSIZE=4 

DELMSH=. 039063 
INCR=1 

ELSEIFCISIZE.EQ.  6)  THEN 
NSIZE=2 
DELMSH=. 0781 
INCR=1 
ELSE 
NSIZE=1 
DELMSH=.  1563 
1NCR=1 
E.NDIF 

WRITEC*,*)' 

WRITEC*,*) 'ANOTHER  INPUT  PARAMETER  IS  THE  SEED  VALUE  OF  THE’ 
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WRITEC*,*) 'RANDOM  NUMBER  GENERATOR.  INPUT  THE  SEED  VALUE  OF  +1.0' 
VRITE(*,*)'OR  -1.0' 

READ(*/OSEED 
WRITE (*,*)' 

WRITE(*,*) 'FINALLY  INPUT  AN  INTEGER  OF  VALUE  1  FOR  TURBULENCE,  OR' 
WRITE(*,*)'0  FOR  NO  TURBULENCE.  ' 

READ(*,*)NYES 
IF(NYES. EQ. 1)  THEN 
WRITE(*,*)' 

WRITE(*,*) ’ INPUT  THE  VALUE  OF  FILTER' 

READ(*,*)FILTER 
END  IF 

WRITEC*,*)' 

WRITEC*,*)' INPUT  THE  VALUE  FOR  CN2' 

READ(*,*)CN2 

WRITEC*,*)' 

c  WRITEC*,*) ' INPUT  THE  RECEIVING  FIELD  SIZE  IN  METERS' 
c  RE ADC*,*) DREG 

WRITEC*,*)'  ' 

WRITEC*,*) ' INPUT  THE  TOTAL  PROPAGATION  DISTANCE  IN  METERS' 

READ;;*,*)Z 

WRITEC*,*)' 

c  WRITEC*,*) ' INPUT  THE  NUMBER  OF  SCREENS  TO  BE  USED  EITHER  FOR' 
c  WRITEC*,*) 'FRESNEL  OR  i'RAUKHOFER  ' 

c  READC*,*)NUMSCR 

IFCSEED. GE. 1. 0)  THEN 
ONECl; 2)='+l’ 

ELSE 

ONECl: 2)=' -1' 

END  IF 

TITLEC20;21)=0.\’EC1:2) 

N2=NR/2 
WVL=. 5E-6 

M=ALOGt:  REALC  .NR ) )  / ALOG C  2 .  ) 

DTRNS=. 3125 
ICNT=0 
C 

C  THE  FOLLOWING  STATEMENT  DETERMINES  WHICH  FORM  OF  FRESNEL  IS  TO  BE 
C  USED 

C 

c  MODE=INTC  C  2*DTRNS*DREC ) / C WVL*Z ) ) 

C 

C  THE  FOLLOWING  SUBROUTINE  CALLED  MCF  DETERMINES  THE  AUTOCORRELATION 

C  OF  THE  APERTURE  WHICH  IS  TO  BE  USED  IN  DETERMINING  THE  ATMOSPHERIC 

C  COHERENCE  LENGTH. 

C 

CALL  MCF  C  FNORM ,  NR ,  M ,  P I ,  NS  I ZE ,  N2 ,  DELMSh :/ 

THIS  SECTION  CREATES  THE  PLANAR  ELECTRIC  FIELD 

DO  40  1=1, NR 
DO  40  J=1..NR 
FIELDRCI,J)=0. 0 
FIELDICI,J)=0. 0 
40  CONTINUE 

DO  41  I=N2-NSIZE+1,N2+NSIZE 
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DO  41  J=N2-NSIZE+1,N2+NSIZE 
FIELDR(I,J)  =  1.  0 
41  CONTINUE 

THIS  SECTION  DETERMINES  THE  SLAB  THICKNESS(ES)  FOR  WHICH  THE 
ELECTRIC  FIELD  IS  PROPAGATED  THROUGH. THIS  IS  VALID  FOR  BOTH 
FORMS  OF  PROPAGATION. 

IF(NR. LE. MODE)  THEN 
DELX=Z/NUMSCR 
ELSE 
DELX=Z 

THE  FOLLOWING  SUBROUTINE  PLACES  A  CURVATURE  ON  THE  WAVEFRONT  TO  BE 
USED  FOR  THE  CONVOLUTION  FORM  OF  FRESNEL  PROPAGATION. 

CALL  QUAD 1 ( NR , P I , DE  LMSH , DX , DY , WVL , DELX ) 

ENDIF 

1000  CONTINUE 

THIS  SECTION  CALLS  OUT  GAUSSIAN  RANDOM  NUMBERS  FOR  THE  REAL  A.ND 
IMAGINARY  PHASE  ARRAYS. 

CALL  GGAUS(.NR,SEED) 

THIS  BEGINS  THE  FILTERING  PROCESS  OF  THE  PHASE  SCREEN 
CALL  FLTR( NR , DELMSH , CN2 , DELX , WVL .FILTER ) 

HERE  THE  INDIRECT  TRANSFORM  IS  BEING  APPLIED  TO  THE  FILTERED 
PHASE  SCREENS. 

CALL  IFTSCR(NR,M, DELMSH) 

THE  FOLLOWING  ARRAYS  USED  IN  THE  FFT  ROUTINES  ARE  ZEROED  OUT  TO 
ENSURE  THAT  UMv’ANTED  VALUES  ARE  NOT  LEFT  IN  THE  ARRAYS. 

DO  990  1=1, NR 
RE(I)=0.  0 
RIM(I)=0.  0 
990  CONTINUE 

THIS  SECTION  DOES  THE  ALGEBRA  NEEDED  TO  MESH  THE  PHASE  SCREEN 
TOGETHER  WITH  THE  ELECTRIC  FIELD 

DO  50  1=1, NR 
DO  50  J=1,NR 

XA=C0S( PHASERC I , J) )*FIELDR( I , J) 

XB=COS( PHASERC I , J) )*FIELDI( I , J) 

XC=SIN( PHASERC I , J) )*FIELDR( I , J) 

XD=S INC PHASERC I , J) )*FIELDI( I , J) 

FIELDRCI,J)=XA-XD 
FIELDICI,J)=XB+XC 
50  CONTINUE 
C 
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HERE  THE  FAST  FOURIER  TRANSFORM  IS  BEING  APPLIED  TO  THE  PERTURBED 
ELECTRIC  FIELD.  FOR  A  DIRECT  TRANSFORM  SIGN=-1.0,  AND  INDIRECT 
TRANSFORM  SIGN=+1. 0. 

SIGN=-1.  0 

CALL  DFTIFT(NR.M,SIGN,DELMSH) 

DO  141  1=1, NR 
DO  141  J=1,NR 

FIELDRC I , J)=FIELDR( I , J)/(NR*DELMSH) 
FIELDI(I,J)=FIELDI(I,J)/(NR*DELMSH) 

141  CONTINUE 


THIS  PORTION  OF  THE  IF  STATEMENT  CORRESPONDS  TO  THE  IMPLEMENTATION 
OF  THE  TRANSFER  FUNCTION  FORM  OF  FRESNEL  PROPAGATION.  THE 
SUBROUTINE  CALLED  TRNSFER  APPLIES  A  QUADRATIC  TO  THE  FIELD. 

IF(NR. LE.  MODE)  THEN 

GALL  TRNSFRC  NR , P I , DX , DY , WVL , DELX , DTRNS ) 

IC.\T=ICNT+1 
SIGN=+1. 0 

CALL  DFTIFT( NR , M , S IGN , DELMSH) 

GO  TO  88S 
ELSE 


THE  SUBROUTINE  CALLED  QUAD2  PUTS  THE  DIFFRACTION  PATTERN  IN  REAL 
SPACE  COORDINATES. 

CALL  QUAD2 ( NR , P I , DX , DY , WVL , DELX) 

ENDIF 

888  CO.NTINUE 


THIS  DO  LOOP  DETERMINES  THE  POWER  SPECTRAL  DENSITY  AND  SETS  IT  UP 
FOR  AN  FFT  TO  DETERMINE  THE  MCF. 

DO  152  1=1, NR 
DO  152  J=1,NR 

FMCF(  I ,  J)=FIELDR(  I ,  J)->*2+FlELDI(  I ,  J)**2 
FILLd  ,J)=0.  0 
152  CO.NTINUE 

ONCE  AGAIN  THE  ARRAYS  ARE  CLEARED  OF  STRAY  VALUES 

DO  910  1=1, NR 
RE(I)=0.  0 
RIM(I)=0.  0 
910  CO.NTI.NUE 

THE  INVERSE  FFT  IS  APPLIED  TO  THE  POWER  SPECTRAL  DENSITY 

SIGN=+1.  0 
DO  901  1=1, NR 
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RE(J)=FMCF(I,J) 

RIM(J)=FILL(I,J) 

911  CONTINUE 

CALL  FFT(M,SIGN,DELMSH) 
DO  921  J=1,NR 
FMCF(I,J)=RE(J) 
FILL(I,J)=RIM(J) 

921  CONTINUE 
901  CONTINUE 

DO  931  J=1,NR 
DO  941  1=1, NR 
RE(I)=FMCF(I,J) 
RIM(I)=FILL(I,J) 

941  CONTINUE 

CALL  FFT(M,SIGN,DELMSH) 
DO  951  1=1, NR 
FMCF(I,J)=RE(I) 
FILL(I,J)=RIM(I) 

951  CONTINUE 
931  CONTINUE 


THIS  SECTION  DETERMINES  THE  MAXIMUM  VALUE  AND  NORMALIZES  THE  MCF 

FLDM=0.  0 
DO  89  1=1, NR 
DO  89  J=1,NR 
XMG=FMCF(I,J) 

IF(XMG.GT. FLDM)  THEN 
FLDM=XMG 
END  IF 

89  CONTINUE 

WRITE(*,*)  FLDM 
PAUSE 

THE  MCF  IS  NORMALIZED  SO  THE  MAX  VALUE  IS  1. 0 

DO  29  1=1, N2 
DO  29  J=1,N2 

FMCF( I,J)=FMCF(I,J)/FLDM 
29  CONTINUE 

THE  ATMOSPHERIC  MCF  IS  DETERMINED  BY  DIVIDING  OUT  THE  APERTURE 
FUNCTION  FROM  THE  PERTURBED  ELECTRIC  FIELD  MCF. 

DO  91  1=1, N2 
DO  91  J=1,N2 

FMCF( I , J)=FMCF( I , J)/FNORM( I , J) 

91  CONTINUE 


DO  87  1=1,1 
DO  87  J=l,20 
WRITE(*,*)FMCF(I,J) 
87  CONTINUE 
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PAUSE 


THIS  SECTION  SETS  UP  THE  ARRAYS  TO  PLOT  THE  2-D  MCF 


NT0T=2*NSIZE 
DO  827  1=1,1 
DO  827  J=1,NT0T 
DISTCJ)=J 
FMAG(J)=FMCF(I, J) 

827  CONTINUE 

DO  83  J=1,NT0T 
WRITE(*,*)FMAG(J) ,DIST( J) 

83  CONTINUE 
PAUSE 

THE  FOLLOWING  SUBROUTINES  ARE  FOR  THE  GRAPHICS  PACKAGE 
THE  ATMOSPHERIC  MCF  IS  BEING  PLOTTED  AT  THIS  POINT 

CALL  VSINIT(18,8.  ,10.  ,0 , 'MCFl.  PLT' , lUNITV, IVID.S) 

CALL  ORIGIN(.5,1.5.0) 

CALL  SGAT,Ei'T^MAG,5.  ,NTOT,n 

CALL  AXIS(0.  ,0.  , 'HCF' ,0,1,1,5.  ,90.  ,FMAG(NTOT+l) ,FMAG(NTOT+2) , .  1 , 1) 
CALL  SCALE(D1ST,5. ,NTOT,l) 

CALL  AXIS(0.  ,0.  ,'  ' ,0,-1, -1,5.  ,0.  ,DIST(NT0T+1) ,DISTCNTOT+2) , . 1,1) 
CALL  LIN£S(DIST,FMAG,NTOT,l,-l,ISYMB,. 1) 

CLOSE(IUNTTV) 

CALL  INT86(ITN0,IREG) 

CALL  MSGCO.  ,0.  ,.  15, 'PRESS  ANY  KEY  TO  CONTINUE*, 0.  ,0,1) 

ANS=GETCf ) 

CALL  GMODECIVID) 


THE  MAGNITUDE  OF  THE  FFT'D  ELECTRIC  FIELD  IS  CALCULATED  IN  ORDER 
TO  PLOT  THE  OUTPUT. 


FMAX=0.  0 
DO  80  1=1, NR 
DO  50  J=1,NR 

FIELDMC I , J)=SQRT(FIELDR(I , J)**2+FIELDI( I , J)**2) 
X=FIELDMCI,J) 

IF(X. GT. FMAX)  THEN 
FMAX=X 
ENDIF 

80  CONTINUE 

THIS  SECTION  BEGINS  THE  CALLING  SEQUENCE  FOR  PLOTTING 

CALL  VSINIT(MON,10.  ,8.  ,0 ,' DITHER. PLT' , lUNITV, IVID,5 ) 
DO  100  1=1, N2 
DO  100  J=1,N2 
IX=J*INCR 
I  Y=  I*  I. NCR 

XMAX=ALOG10(FIELDM( I ,J)/FMAX) 


onnn  onoo 


INDEX=8*ABS(XMAX)+1 
IF( INDEX. GE. 21)  THEN 
ICOLOR=0 
ELSE 

ICOLOR=NDEX( INDEX) 

END  IF 

CALL  PIXEL(IX,IY,ICOLOR) 

100  CONTINUE 

CALL  MSG(0.  ,1-  ,•  15, TITLE, 0.  ,0,0) 

CLOSE (I UNI TV) 

C  CALL  INT86(1TN0,IREG) 

CALL  MSG(0. ,0. ,. 15, 'PRESS  ANY  KEY  TO  CONTINUE ',0. ,0,0) 

ANS=GETC( ) 

CALL  GMODE(IVID) 

THIS  IF  STATEMENT  QUES  THE  PROGRAM  TO  START  OVER  AGAIN  IF 
FRAUNHOFER  OR  THE  CONVOLUTION  FORM  OF  FRESNEL  ARE  USED 

IF(NR. GT. MODE)  THEN 
GO  TO  999 
ENDIF 

THIS  IF  STATEMENT  QUES  THE  PROGRAM  TO  FULLY  COMPLETE  PROPAGATION 
THROUGH  ALL  THE  SLABS  IN  THE  TRANSFER  FUNCTION  FORM  OF  FRESENL 

IF( ICNT. NE.  NUMSCR)  THEN 
GO  TO  1000 
ENDIF 
GO  TO  999 
E.ND 

SUBROUTINE  QUAD1(NR.PI ,DELMSH,DX,DY,WVL,DELX) 

COMMON  /BLK2/  FIELDRC 256 , 256) ,F1ELDI( 256 , 256) 

D.X=DELMSH 
DY=DELMSK 
MID=(NR/2)+l 
DO  273  1=1, NR 

X=(I-(MID-"2-l))*DX 
DO  273  J=1,NR 

Y=( J-(MID*2-1))*DY 
THETA=P  I  ( (  (  X*X )  +  (.  Y* Y  )  )  /  (  WVL*DELX) ) 

XX=FIELDR( I , J)*COS(THETA) 

YY=FIELDR( I , J)*SIN(THETA) 

ZZ=FIELDI(I,J)*COS(THETA) 

W=F  I E  LD I  ( I ,  J )  *  S I N  ( THET A ) 

FIELDRC  I,  J)=XX-W 
FIELDI(I,J)=YY+ZZ 
273  CONTINUE 
RETURN 
END 

SUBROUTINE  QUAD2(NR,PI ,DX,DY,WVL,DELX) 

COMMON  /BLK2/  FIELDR(256,256) ,FIELDI(256,256) 

DX2=DX-"-WVL-'--DELX 
DY2=DY*v,'VL*DELX 
MID=NR/2+l 


S') 
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DO  274  1=1, NR 
Y=(I-MID)*DY2 
DO  274  J=1,N’R 
X=(J-MID)*DX2 

PHI=PI*(  ( (X*XH(y*Y)  )/(WVL*DELX) ) 
CX=FIELDR( I , J)*COS( PHI ) 
CY=FIELDR(I,J)*SIN(PHI) 
CZ=FIELDI(I,J)*COS(PHI) 
CW=FIELDI(I,J)*SIN(PHI) 

CBR=CX-CW 

CBI=CY+C7 

F I ELDR ( I , J ) = ( C  B I / ( W VL*DE  LX ) ) 
FIELDI(I,J)=-1.*CCBR/(WVL*DELX)) 

274  CONTINUE 
RETURN 
END 
C 

SUBROUTINE  TRNSFR(NR,PI ,DX,DY,WVL,DELX,DTRNS) 
COMMON  /BLK2/  FIELDR( 256 , 256) ,FIELDIC256 , 256) 
DX=DTRNS/NR 
Dy=DTRNS/NR 
MID=NR/2+l 
DO  275  1=1, NR 
FY=(I-MID)*DY 
DO  275  J=1,NR 
FX=(J-MID)*DX 

FEE=- 1.  *PI--'WVL*DELX*(  (FX*FX)+(FY*FY) ) 

GX=F IE  LDR ( I , J ) *COS ( FEE ) 

GY=i  lELDRC  I ,  J)->SIN(FEE) 
GZ=FIELDI(I,J)*COS(FEE) 

GW=FIELDI(I ,J)*SIN(FEE) 

FIELDR(I,J)=GX-GW 
FIELDI(I,J)=GY+GZ 
273  CONTINUE 
RETURN 
E.'vD 


SUBROUTINE  FLTR( NR , DELMSH , CN2 , DELX , WVL , FI LTER) 
C  NOTE:  DELMSH^'-NR  IS  THE  LARGEST  APERTURE  SIZE 
COMMON  /BLK3/  PHASER(256,256) ,PHASEI(256,256) 
PI=3. 141592653589792 

pov;er=-ii.  /6. 

TPI=2. *PI 
N2=NR/2 
NPIVOT=N2+l 
LAST=NPIVOT+l 

DLKAPA=( TPI/ ( NR*DELMSH) )**POWER 
FACTOR=SQRT( (TPI**3)*.  033*CN2*DELX/(WVL**2) ) 
FUDGE=DLKAPA*FACTOR 
DO  100  I=l,NPIVOT 
EYE=REAL( I ) 

EYE2=EYE-”EYE 
DO  100  J=1,NR 
IF(J. LE. NPIVOT)  THEN 
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WHY=REAL(J) 

ELSE 

VHY=REAL(NR-J+2) 

END  IF 

XKAPPA=SQRT( EYE  2+WHY*WHY ) 

AKAPPA=XKAPPA**FILTER 
PHASER( I , J)=PHASER( 1 , J)*FUDGE*AKAPPA 
PHASE I ( I , J ) =PHASE I ( I , J ) *FUDGE*AKAPPA 
100  CONTINUE 

DO  110  I=LAST,  NR 
EYE=REAL(NR-I+2) 

EYE2=EYE*EYE 
DO  110  J=1,NR 
IF(J. LE.NPIVOT)  THEN 
WHY=REAL(J) 

ELSE 

WHY=REAL(NR-J+2) 

ENDIF 

XKAPPA=SQRT( EYE2+WHY*WHY) 

AKAPPA=XKAPPA**F I LTER 
PHASERC I , J)=PHASER( I , J)*FUDCE*AKAPPA 
PHASEK I , J)=PHASEI( I , J)*FUDGE*AKAPPA 
110  CONTINUE 
C 

PHASERC 1,1)=0. 0 
PHASEK  1,1)=0.  0 

NOTE:  .  03441=TPI**-11. /6.  TO  TRANSFORM  FROM  KAPPA  TO  FREQ  SPACE. 
DO  11  1=1, NR 
DO  11  J=1,.\R 

FHASER(I,J)=PHASER(I,J)*. 03441 
PKASEi(I,J)=PHASEI(I,J)*. 03441 
11  CO.NTI.NUE 
RETURN 
END 

SUBROUTINE  IFTSCR(NR,M,DELMSH) 

COMMON  /BLKl/  REC 256) ,RIM(256) 

COM.MON  /BLK3/  PHASERC  256 ,256)  ,PHASEI( 256, 256) 

PI=3. 141592653589792 
TPI=2. *PI 
SIG.\=-1.  0 
DO  20  1=1, NR 
DO  21  J=1,NR 

RECJ)=PHASERCI,J) 

RIMCJ)=PHASE1CI,J) 

21  CONTINUE 

CALL  FFTCM,SIGN,DELMSH) 

DO  22  J=1,NR 

PHASERC I, J)=RECJ) 

PHASEICI,J)=RIMCJ) 

22  CONTINUE 
20  CONTINUE 

DO  30  J=1,NR 
DO  31  1=1, NR 

RECI)=PHASERCI,J) 
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RIM(I)=PHASEI(I,J) 

31  CONTINUE 

CALL  FFT(M,SIGN,DELMSH) 

DO  32  1=1, NR 

PHASER(I,J)=RE(I) 

PHASEI(I,J)=RIM(I) 

32  CONTINUE 
30  CONTINUE 

RETURN 

END 

SUBROUTINE  DFTIFT( NR , M , S I GN , DELMSH ) 

COMMON  /BLKl/  RE(256) ,RIM(256) 

COMMON  /BLK2/  FIELDR(256,256) ,FIELDI(256,256) 
DO  60  1=1, NR 
DO  61  J=1,NR 

RE(J)=FIELDR(I,J) 

RIM(J)=FIELDI(I,J) 

61  CONTINUE 

CALL  FFT(M, SIGN, DELMSH) 

DO  62  J=1,NR 

FIELDR(I,J)=RE(J) 

FIELDI(I,J)=RIM(J) 

62  CONTINUE 
60  CONTINUE 

DO  70  J=1.NR 
DO  71  1=1, NR 

RE(1)=FIELDR(I,J) 

RIM(I)=FIELDI(I,J) 

71  CONTINUE 

CALL  FFT(H, SIGN, DELMSH) 

DO  72  1=1. NR 

FTELD?.ri,J)=RE(I) 

FIELDI(I,J;=RIM(I) 

72  CONTINUE 
70  CONTINUE 

RETURN 

END 


SUBROUTINE  FFT(M, SIGN, DELMSH) 

COMMON  /ELKl/  RE( 256) ,RIM( 256) 

PI=3.  141592653589792-'’^SIGN 

N=2**M 

N1=N-1 

J=1 

DO  200  1=1, N1 
IFCI.LT.  J)  THEN 
T=RE(J) 

RE(J)=RE(I) 

RE(I)=T 

T=RIM(J) 

RIMCJ)=RIM(I) 

RIM(I)=T 
END  IF 
K=N/2 
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DO  201  WHILECK. LT. J) 

J=J-K 

K=K/2 

201  CONTINUE 
J=J+K 

200  CONTINUE 
LE=1 

DO  202  L=1,M 
LE1=LE 
LE=LE+LE 
URE=1. 

UIM=0. 

ANG=PI/LE1 

VRE=COS(ANG) 

WIM=SIN(ANG) 

DO  203  J=1,LE1 
DO  204  I=J,N,LE 
IP=I+LE1 

TRE=RE( IP)*URE-RIM( IP)*UIM 

TIM=RE( IP)*UiM+RIMC IP)*URE 

RE(IP)=RE(n-TRE 

RIM(IP)=RIM(I)-TIM 

RE(I)=RE(I)+TRE 

RIM(I)=RIM(I)+TIM 

204  CONTINUE 
T=URE--nvRE-UIN*VIM 

U I  M=URE’nv'I  M+U I WRE 
URE=T 

203  CONTINUE 

202  CONTINUE 
IFCSIGN.GT. 0. 0)  THEN 
PTS=1.  0/(N--'-DELMSH) 

BO  205  1=1, N 

RE(I)=RE(I)^-PTS 

RIM(I)=R1M(I)*PTS 

205  CONTINUE 
ENPIF 
RETURN 
END 

C 

C 

SUBROUTINE  GGAUS(NR,SEED) 

COMMON  /ELK3/  PHASER(256,256) ,PHASEI(256,256) 
DO  300  1=1, NR 
DO  300  J=1,NR 
301  Vl=2. *RAN(SEED)-1 
V2=2.  *RAN(SEED)-1 
S=V1*V1+V2^>V2 
IF(S. GE. 1. 0)  GO  TO  301 
SCALE=SQRT(-2.  *ALOG(S)/S) 

X1=V1-“'SCALE 
C  X2=V2*SCALE 

PHASER(I,J)=X1 
PHASEI(I,J)=0.  0 
300  CONTINUE 
RETURN 
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END 

C 

SUBROUTINE  MCF( FNORM , NR , M , P I , NS I ZE , N2 , DELMSH ) 
COMMON  /BLKl/  RE(256) ,RIM(256) 

COMMON  /BLK2/  FIELDR( 256 , 256) ,FIELDI( 256 , 256) 
DIMENSION  FN0RMC256,256) 

THIS  SECTION  CREATES  THE  PLANAR  ELECTRIC  FIELD 

DO  39  1=1, NR 
DO  39  J=1,NR 
FIELDR(I,J)=0.  0 
FIELDI(I,J)=0.  0 
39  CONTINUE 

DO  45  I=N2-NSI2E+1,N2+NSI2E 
DO  45  J=N2-NSIZE+1,N2+NSIZE 
FIELDR(I,J)  =  1.  0 
45  CONTINUE 


SIGN=-1.  0 

CALL  DFT I FT( NR , M , S I GN , DE LMSH ) 


DO  80  1=1, NR 
DO  80  J=1,NR 

FN0RM( I , J)=FIELDR( I , J)**2+FIELDI( I , J)**2 
80  CONTINUE 

Du  120  1=1, NR 
DO  123  J=1,NR 
FIELDI(I,J)=0.  0 
123  CONTINUE 

SIGN=+1.  0 

DO  90  1=1, NR 
DO  91  J=1,NR 
RE(J)=FNORM(I,J) 

RIM(J)=FIELDI(I,J) 

91  CONTINUE 

CALL  FFT(M, SIGN, DELMSH) 

DO  92  J=1,NR 
FNORM(I,J)=RE(J) 

FIELDI(I,J)=RIM(J) 

92  CONTINUE 
90  CONTINUE 

DO  93  J=1,NR 
DO  94  1=1, NR 
RE(I)=FNORM(I,J) 

RIM(I)=FIELDI(I,J) 

94  CONTINUE 

CALL  FFT(M, SIGN, DELMSH) 

DO  95  1=1, NR 
FNORM(I,J)=RE(I) 

FIELDI(I,j)=RIM(I) 


o  o 


I 


95  CONTINUE 
93  CONTINUE 


FLDM=0.  0 
DO  88  1=1, NR 
DO  88  J=1,NR 
XMG=FNORM(I,J) 

IF(XMG. GT. FLDM)  THEN 
FLDM=XMG 
END  IF 

88  CONTINUE 

DO  89  1=1, N2 
DO  89  J=1,N2 

FNORMC 1 , J)=FNORM( I , J) /FLDM 

89  CONTINUE 

RETURN 

END 
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